library(DEP)
library(DEqMS)
library(ggpubr)
library(tidySummarizedExperiment)
library(tidyverse)

kegg_mva <-
  c('Acat1','Acat2','Hmgcs1','Hmgcs2','Hmgcr','Mvk','Pmvk','Mvd','Idi1','Idi2','Fdps') |>
  str_to_upper()

meta_xy <- read_delim('mission/FPP/xiangya_sle_ms/ms_sample_meta.txt')

meta_wst <- meta_xy |>
  mutate(label = str_c(batch, '_', label),
         condition = case_when(str_detect(Group, 'ann') ~ 'aSCLE',
                               str_detect(Group, 'papu') ~ 'pSCLE',
                               .default = Group)) |>
  select(label, condition) |>
  group_by(condition) |>
  mutate(replicate = seq_along(label)) |>
  ungroup()

wst_design <- colnames(wstlk)[3:47] |>
  as_tibble() |>
  rename(label = value) |>
  left_join(meta_wst) |>
  mutate(batch = str_extract(label, 'F.'))

# westlake omics mtx --------
wstlk <- read_csv('mission/FPP/xiangya_sle_ms/xiangya_matrix.csv') |>
  separate_wider_delim(...1, delim = '_',
                       names = c('uniprot','protein'),
                       too_few = 'align_start', too_many = 'drop')

## uniprot suffix is not always the same as symbol
## uniprot:symbol pair will return 1:many
wstlk_symbol <- wstlk$uniprot |>
  clusterProfiler::bitr(fromType = 'UNIPROT',
                        toType = 'SYMBOL',
                        OrgDb = 'org.Hs.eg.db') |>
  distinct(UNIPROT, .keep_all = TRUE) |>
  rename(uniprot = UNIPROT) |>
  right_join(wstlk) |>
  as_tibble()

wstlk_symbol |>
  filter(is.na(SYMBOL) & protein == 'NA')

wstlk |>
  pivot_longer(3:last_col()) |>
  ggplot(aes(name, value)) +
  geom_boxplot() +
  coord_flip()
  
se_wstlk <- wstlk_symbol |>
  make_unique('SYMBOL','uniprot') |>
  make_se(4:48, wst_design)

se_wstlk

## remove 5 outlier --------
se_wstlk <- se_wstlk |>
  filter(!str_detect(label, '133C|3_132C|2_132N'))

se_wstlk |> plot_frequency()
se_wstlk |> plot_numbers()

se_wstlk |> plot_pca(indicate = 'batch', n=6990)

# Filter for proteins that are identified in all replicates of at least one condition
# 6990 -> 5358
se_nomiss <- se_wstlk |> filter_missval()

se_nomiss |> plot_frequency()
se_nomiss |> plot_numbers()
se_nomiss |> plot_pca(indicate = 'batch', n = 5000)

## normalize ---------
se_norm <- se_nomiss |> normalize_vsn()

se_nomiss |> plot_normalization(se_norm)
se_norm |> plot_pca(indicate = 'batch', n = 5000)

## impute missing ---------
## se_nomiss already have no missing value
se_wstlk |> plot_detect()

se_miss <- se_wstlk |> normalize_vsn()
se_miss |> plot_pca(indicate = 'batch', n = 5000)

# for Missing At Random use knn impute
se_imp <- se_miss |>
  impute(fun = "min")

se_miss |>
  plot_imputation(se_imp)

## test diff ------
se_diff <- se_imp |>
  test_diff(control = 'NC') |>
  add_rejections(lfc = log2(1.25))

se_diff |> plot_pca(indicate = 'batch')

### Pearson cor matrix -------
se_diff |> plot_cor(pal = "Reds")

se_diff |> plot_heatmap()

# strange parsing warning & wrong plot?
se_diff |> plot_single('HMGCR')

se_diff |>
  plot_volcano(contrast = "SLE_vs_NC")

res_diff <- se_diff |>
  get_results() |>
  as_tibble()

res_diff |>
  filter(name == 'HMGCR')

## test diff (no miss) ----------
se_nomiff <- se_norm |>
  test_diff(control = 'NC') |>
  add_rejections(lfc = log2(1.25))

se_nomiff@colData <- se_nomiff@colData |>
  as.data.frame() |>
  mutate(batch = str_extract(label, 'F.'))

se_nomiff |> plot_pca(indicate = 'batch')

se_nomiff |> plot_cor(pal = "Reds")

se_nomiff |> plot_heatmap()

se_nomiff |> plot_single('HMGCR')

se_nomiff |>
  plot_volcano(contrast = "SLE_vs_NC")

res_nomiff <- se_nomiff |>
  get_results() |>
  as_tibble()

res_nomiff |>
  filter(SLE_vs_NC_p.adj < .05)

## welch t test -----------
wst_mva <- wstlk |>
  #filter(protein %in% kegg_mva) |>
  select(-uniprot) |>
  pivot_longer(-1, names_to = 'label') |>
  left_join(wst_design) |>
  select(-replicate)

wst_welch <- 
wst_mva |>
  filter(condition %in% c('NC','SLE') & !is.na(value) & protein != 'NA') |>
  filter(!str_detect(label, '133C|3_132C|2_132N')) |>
  reframe(value = list(value), .by = c(protein, condition)) |>
  pivot_wider(names_from = condition, values_from = value) |>
  rowwise() |>
  filter(length(SLE) >2 & length(NC) > 2) |>
  mutate(pval = t.test(SLE, NC)$p.value,
         log2fc = log2(t.test(SLE, NC)$estimate[1] / t.test(SLE, NC)$estimate[2])) |>
  select(protein, pval, log2fc) |>
  ungroup() |>
  mutate(padj = p.adjust(pval, method = 'fdr'))

wst_welch |>
  ggplot(aes(log2fc, -log10(padj))) +
  geom_point()

wst_welch |>
  filter(log2fc < -log2(1.25) & padj < .05)

## DEqMS -------------
# define group meta
wst_fct <- wst_design$condition |>
  fct()

wst.model.mat <- model.matrix(~ 0 + wst_fct)

colnames(wst.model.mat) <- wst.model.mat |>
  colnames() |>
  str_remove('wst_fct')

# define contrast to test
wst.contrast <- makeContrasts('SLE-NC', levels = wst.model.mat)

wstlk.log <- wstlk |>
  na.omit() |>
  select(-protein) |>
  mutate(across(-1, log2)) |>
  column_to_rownames('uniprot')

wstlk.limma <- wstlk.log |>
  lmFit(wst.model.mat) |>
  contrasts.fit(wst.contrast) |>
  eBayes()

head(wstlk.limma$coefficients)

# without PSM correction, can only use limma::topTable to get res
limma_res <- wstlk.limma |>
  topTable(coef = 1, n = Inf) |>
  as_tibble()

limma_res |> filter(adj.P.Val < .05)

# openms res-----------
# msgf
msgf <- read_tsv('~/append-ssd/nextflowing/quantms_xy_b123_msgf_canon/proteinquantifier/xiangya.b123all.sdrf_openms_design_protein_openms.csv', comment = '#')

msgf <- msgf |>
  mutate(protein = str_remove(protein, '.+/')) |>
  separate_wider_delim(protein, '|',
                       names = c('foo','uniprot'), too_many = 'drop') |>
  select(-c(foo, n_proteins, protein_score, contains('ratio')))

msgf <- msgf$uniprot |>
  clusterProfiler::bitr(fromType = 'UNIPROT',
                        toType = 'SYMBOL',
                        OrgDb = 'org.Hs.eg.db') |>
  rename(uniprot = UNIPROT) |>
  left_join(msgf) |>
  as_tibble() 

opms_design <- meta_xy |>
  filter(batch == 'F1') |>
  mutate(label = openms.id |> as.character(),
         condition = Group,
         .keep = 'none') |>
  group_by(condition) |>
  mutate(replicate = seq_along(condition))

opms_design

opb1_se <- b1e2 |>
  make_unique('SYMBOL','uniprot') |>
  make_se(3:17, opms_design)

opb1_se |> plot_frequency()

opb1_se |> plot_numbers()

# Filter for proteins that are identified in all replicates of at least one condition
# 7110 -> 7110
opb1_se_m1 <- opb1_se |>
  filter_missval()

opb1_se_m1

## normalize ---------
opb1_se_norm <- opb1_se |> normalize_vsn()

opb1_se |> plot_normalization(opb1_se_norm)

## impute missing ---------
## openms result seem like Missing Not At Random (MNAR)
opb1_se |> plot_detect()

opb1_se_minprob <- opb1_se_norm |>
  DEP::impute(fun = "MinProb")

opb1_se_norm |>
  plot_imputation(opb1_se_minprob)

opb1_se_man <- opb1_se_norm |>
  DEP::impute(fun = "man")

opb1_se_norm |>
  plot_imputation(opb1_se_man)

opb1_se_zero <- opb1_se_norm |>
  DEP::impute(fun = "zero")

opb1_se_zero |>
  plot_imputation(opb1_se_minprob)

opb1_se_qrilc <- opb1_se_norm |>
  DEP::impute(fun = "QRILC")

opb1_se_qrilc |>
  plot_imputation(opb1_se_minprob)

## test diff ------
opb1_se_minprob <- opb1_se_minprob |>
  test_diff(control = 'NC') |>
  add_rejections()

opb1_se_man <- opb1_se_man |>
  test_diff(control = 'NC') |>
  add_rejections()

opb1_se_zero <- opb1_se_zero |>
  test_diff(control = 'NC') |>
  add_rejections()

opb1_se_qrilc <- opb1_se_qrilc |>
  test_diff(control = 'NC') |>
  add_rejections()

## Pearson cor matrix -------
opb1_se_dep |> plot_cor(pal = "Reds")

opb1_se_minprob |>
  plot_volcano(contrast = "SLE_vs_NC", adjusted = T)

opb1_se_man |>
  plot_volcano(contrast = "SLE_vs_NC", adjusted = T)

opb1_se_zero |>
  plot_volcano(contrast = "SLE_vs_NC", adjusted = T)

opb1_se_qrilc |>
  plot_volcano(contrast = 'SLE_vs_NC', adjusted = T)

opb1_se_zero |> plot_single('MVK')

opb1_se_minprob |> get_results() |>
  as_tibble() |>
  dplyr::count(contains('signific'))

# from msstats -----------
b13ms_prot_mtx <- read_csv('mission/FPP/xiangya_sle_ms/sle_ms_batch13_e2_wide_mtx.csv')

b13col <- b13ms_prot_mtx |>
  colnames()

b13design <- tibble(label = b13col[3:32]) |>
  mutate(condition = str_extract(label, 'NC|.+LE')) |>
  group_by(condition) |>
  mutate(replicate = seq_along(label))

b13dep <- b13ms_prot_mtx |>
  make_unique('SYMBOL','UNIPROT') |>
  make_se(3:32, b13design)

b13dep |> plot_frequency()

# default filter for prot that are identified in all replicates of at least one condition
b13filt <- b13dep |>
  filter_missval()

b13filt |> plot_frequency()

b13dep |> plot_numbers()

## remove outliers -------
meta_xy |>
  filter(label == '133C' | 
           batch == 'F3' & label == '132C' |
           batch == 'F2' & label == '132N')

b13dep <- b13dep |>
  filter(!(.sample %in% c('SLE_4','DLE_4','annularSCLE_1')))

b13dep |> plot_numbers()

b13dep |> plot_frequency()

b13filt <- b13dep |>
  filter_missval()

b13filt |> plot_frequency()

# 2nd impute

# test diff on already normalized & imputed data
# default mode is `control`
b13diff <- b13filt |>
  test_diff(control = 'NC')

# default p<0.05 lfc>1
b13res <- b13diff |>
  add_rejections()

b13_msstats_dep <- b13res |>
  get_results() |>
  as_tibble()

b13_msstats_dep |>
  filter(SLE_vs_NC_p.adj < .05) |>
  slice_min(SLE_vs_NC_ratio) |>
  select(SLE_vs_NC)

b13res |> plot_volcano(contrast = 'SLE_vs_NC')

b13_slevshc <- b13res |>
  rowData() |>
  as_tibble() |>
  select(c(SYMBOL, UNIPROT, contains('SLE_vs')))

b13res |> plot_pca()

b13res |> plot_heatmap()

b13res |> plot_single('HMGCR')
